R Notebook

library(StereoMorph)
library(tidyverse)

Loading tidyverse: ggplot2
Loading tidyverse: tibble
Loading tidyverse: tidyr
Loading tidyverse: readr
Loading tidyverse: purrr
Loading tidyverse: dplyr
Conflicts with tidy packages ----------------------------------------------------
filter(): dplyr, stats
lag():    dplyr, stats

library(geomorph)

Loading required package: rgl
Loading required package: ape

library(listviewer)
library(stringr)

(localities<-read_csv("Marks_Mox - Sheet1.csv"))

Parsed with column specification:
cols(
  CatalogNumber = col_integer(),
  ScientificName = col_character(),
  `# of individuals` = col_character(),
  Locality = col_character(),
  Country = col_character(),
  YearCollected = col_integer(),
  Photo = col_character(),
  species = col_character(),
  drainage = col_character(),
  river = col_character()
)

localities %>%
  select(CatalogNumber,species,drainage,river)%>%
  rename(cat_num=CatalogNumber)->loc

#files<-paste0("Mox_images/shapes/TU198404_",1:6,"_L.txt")
a<-readShapes(file = "Mox_images/shapes/",fields=c("landmarks.scaled","curves.scaled"))

Put landmark data into a tibble.

data<-tibble(num=attr(a$landmarks.scaled,"dimnames")[[3]],
       fixedlm=array_branch(a$landmarks.scaled,margin = 3),
       c_body_ant=map(a$curves.scaled,"body_ant"),
       c_body_post=map(a$curves.scaled,"body_post"),
       c_opercle=map(a$curves.scaled,"opercle"))

data %>%
  mutate(cat_num=str_extract(num,"[0-9]+"))->data

data$cat_num<-as.integer(data$cat_num)
(data<-left_join(data,loc))

Joining, by = "cat_num"

convert curves into landmarks, evenly spaced along curves. Bind fixed and semi-landmarks together and remove duplicates.

data %>%
  mutate(body_ant=map(c_body_ant,~pointsAtEvenSpacing(.x,n=10)))%>%
  mutate(body_post=map(c_body_post,~pointsAtEvenSpacing(.x,n=10)))%>%
  mutate(opercle=map(c_opercle,~pointsAtEvenSpacing(.x,n=5))) %>%
  mutate(land_marks=pmap(list(fixedlm,body_ant,body_post,opercle),rbind))%>%
  mutate(land_marks=map(land_marks,~unique(.x)))->data

Convert list (and bind several arrays together) using sapply()

new_a<-sapply(data$land_marks, I, simplify="array")

Generate plot to aid in defining sliding, semi-landmarks using AUTO mode of define.sliders(). There are 20 fixed landmarks, curves are found between landmarks 1 and 2 (anterior, dorsal body), 3 and 4 (posterior, dorsal body), and 10 and 11 (opercle)

dd<-as.data.frame(new_a[,,1])
dd$label<-1:length(dd$V1)
ggplot(dd,aes(V1,V2))+
  geom_point(alpha=0.7)+
  geom_text(label=dd$label,check_overlap = F,nudge_x = 1,size=3)

Generate semi-landmarks matrix for gpagen using define.sliders().

curves<-rbind(define.sliders(c(1,29:36,2)),
              define.sliders(c(3,21:28,4)),
              define.sliders(c(10,37:39,11)))

gpa<-gpagen(new_a,curves = curves)

  |                                                                             
  |                                                                       |   0%
  |                                                                             
  |==============                                                         |  20%
  |                                                                             
  |============================                                           |  40%
  |                                                                             
  |===========================================                            |  60%
  |                                                                             
  |=========================================================              |  80%
  |                                                                             
  |=======================================================================| 100%

gpa

Call:
gpagen(A = new_a, curves = curves) 



Generalized Procrustes Analysis
with Partial Procrustes Superimposition

20 fixed landmarks
19 semilandmarks (sliders)
2-dimensional landmarks
6 GPA iterations to converge
Minimized squared Procrustes Distance used


Consensus (mean) Configuration

                 X             Y
 [1,] -0.189292837 -0.0188987033
 [2,]  0.029581528  0.0433555969
 [3,]  0.103513405  0.0240097538
 [4,]  0.291057191  0.0149602944
 [5,]  0.292731582 -0.0399869752
 [6,]  0.222333149 -0.0435116367
 [7,]  0.181608223 -0.0589950374
 [8,]  0.055101051 -0.0671944886
 [9,] -0.083371192 -0.0516103724
[10,] -0.094910572 -0.0507079052
[11,] -0.101084694  0.0135698686
[12,] -0.147023141 -0.0044754231
[13,] -0.157494327 -0.0047139364
[14,] -0.154682252  0.0023113931
[15,] -0.146676471  0.0066572908
[16,] -0.138713301  0.0041156807
[17,] -0.136003081 -0.0041973511
[18,] -0.138979090 -0.0116328834
[19,] -0.147209135 -0.0150931918
[20,] -0.154832756 -0.0115755409
[21,] -0.172234467  0.0001636571
[22,] -0.150291550  0.0134418372
[23,] -0.126774867  0.0226499323
[24,] -0.102904749  0.0295011512
[25,] -0.077958549  0.0350546513
[26,] -0.051693480  0.0387790044
[27,] -0.024679231  0.0411041028
[28,]  0.003760085  0.0422855421
[29,]  0.124162947  0.0208084767
[30,]  0.145958581  0.0178719752
[31,]  0.167119641  0.0151792999
[32,]  0.188043939  0.0127899438
[33,]  0.208877947  0.0106387953
[34,]  0.229624137  0.0089623537
[35,]  0.250128884  0.0083247394
[36,]  0.270922320  0.0107404774
[37,] -0.087803161 -0.0356431614
[38,] -0.087333514 -0.0180416327
[39,] -0.092578192 -0.0009975782

plotAllSpecimens(gpa$coords,mean=F)

plot PCA

PCA <- plotTangentSpace(gpa$coords,warpgrids = F)

ggplot(pca, aes(PC1, PC2, color = id2)) +
  geom_polygon(data = hull,aes(PC1,PC2,group=id,fill=id),alpha=0.3) +
  geom_point()+
  facet_wrap(~id2)

M<-mshape(gpa$coords)
PC<-PCA$pc.scores[,1:2]
preds<-shape.predictor(gpa$coords,x=PC,Intercept = FALSE,
                       pred1=c(-0.05,0.04)) 
                       
GP<-gridPar(pt.size=0.5,tar.pt.size=0.5,n.col.cell = 50)      
plotRefToTarget(M,preds$pred1,mag = 2,method = "vector",gridPars = GP)